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We explore the properties of bosonic atoms loaded into the d-bands of an isotropic square optical 
lattice. Following the recent experimental success reported in [Y. Zhai et al., Phys. Rev. A 87, 
063638 (2013)], in which populating d bands with a 99% fidelity was demonstrated, we present 
a theoretical study of the possible phases that can appear in this system. Using the Gutzwiller 
ansatz for the three d band orbitals we map the boundaries of the Mott insulating phases. For not 
too large occupation, two of the orbitals are predominantly occupied, while the third, of a slightly 
higher energy, remains almost unpopulated. In this regime, in the superfluid phase we find the 
formation of a vortex lattice, where the vortices come in vortex/anti-vortex pairs with two pairs 
locked to every site. Due to the orientation of the vortices time-reversal symmetry is spontaneously 
broken. This state also breaks a discrete Z 2 -symmetry. We further derive an effective spin-1/2 
model that describe the relevant physics of the lowest Mott-phase with unit filling. We argue that 
the corresponding two dimensional phase diagram should be rich with several different phases. We 
also explain how to generate antisymmetric spin interactions that can give rise to novel effects like 
spin canting. 

PACS numbers: 45.50.Pq, 03.65.Vf, 31.50Gh 


I. INTRODUCTION 

Interest in systems of cold atoms in optical lattices 
has greatly increased during the last decade [1] partly 
because of their versatility as simulators of quantum 
systems. More precisely, the flexibility, control, and 
cleanness of these systems have led to numerous realisa¬ 
tions of paradigm many-body models of condensed mat¬ 
ter physics. The field was greatly stimulated with ex¬ 
perimental explorations of the Mott insulator-superfluid 
transition of the Bose-Hubbard (BH) model [2] following 
the theoretical proposal of Ref. [3]. Today the list of 
achievements is long, but to mention just a few: single¬ 
site resolved detection [3], simulation of magnetism and 
spin models [5] , realising synthetic magnetic fields [S] and 
topological states [7] , demonstration of a fermion Mott in¬ 
sulator [5] , as well as various dynamical studies like equi¬ 
libration and Lieb-Robinson spread of correlations [S]. 
Common to all the above examples is that the essen¬ 
tial physics appear on the lowest energy band, the s 
band. In particular, for spinless particles on the s band, 
the onsite atomic states are not degenerate. We know, 
however, that many interesting phenomena in condensed 
matter theories have their origin in the degeneracies of 
onsite states/orbitals. For example, understanding ’or¬ 
bital physics ’ m is essential for giving a full descrip¬ 
tion of superfluid properties of ^He m and the metal- 
insulating transitions in metal oxides [12) . Spurred by 
this, in recent years the first steps towards controlled 
studies of orbital physics with cold atoms have been 
taken [13]. 

Degenerate orbital atomic states appear naturally on 
excited bands of optical lattices, and it was predicted that 


bosonic atoms loaded into the first excited, the p bands, 
of a two dimensional (2D) square optical lattice give rise 
to interesting physics beyond that found on the lowest 
band M- In particular, the superfluid order parameter 
at zero temperature is complex valued and time-reversal 
symmetry is spontaneously broken |14Hltij . Also the in¬ 
sulating phases provide rich physics with the possibili¬ 
ties to study exotic quantum magnetic phases ininH]. 
In different lattice geometries, the system may display 
’stripped phases ’ HH, and in 3D the state of bosons oc¬ 
cupying an isotropic cubic lattice become frustrated both 
in the insulating as well as in the superfluid phase [201121) . 
Fermionic systems in the p band are also known to fea¬ 
ture rich physics, with a plethora of different phases [22) . 

Experimentally, populating the p bands with bosonic 
atoms was first realised by accelerating the lattice po¬ 
tential in such a way that the system traversed through 
a Landau-Zener transition non-adiabatically [23) . How¬ 
ever, a more thorough experimental analysis of the co¬ 
herence and life-time of p-band bosons was performed by 
Muller et al. [M] using a two-laser Raman loading tech¬ 
nique from an insulating s-band state. In particular, it 
was found that the atoms on the p band relax very rapidly 
(on the order of a few tunneling times), and that the 
life-time is relatively long in comparison to the tunnel¬ 
ing time. This enables coherence to be established from 
the initial Mott insulator and the possibility to detect 
dynamical processes. More recently, a 2D superlattice 
was used in order to prepare bosonic atoms in hybridised 
states composed of s and p orbitals [5S]. In this exper¬ 
iment, and in a follow up one [26) . clear evidences of a 
complex order parameter were presented. The same load¬ 
ing method was further benchmarked by loading atoms 
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into /-orbital states i.e. the third set of excited 
bands. 

A completely different approach to initialise bosonic 
atoms in excited bands was proposed and demonstrated 
by Zhai et al. [5S]. There the atoms start in a zero 
momentum eigenstate (held in a broad harmonic trap) 
and then a sequence of ‘on/off pulses’ are applied to the 
atoms, where the ‘on pulse’ consists of a standing wave 
which couple the zero momentum states to higher mo¬ 
mentum states by photon recoils. With this method the 
atoms were loaded into the d band with a fidelity as high 
as 99%, which is much better than what was reported for 
the p bands (around 80%). This technique is, however, 
more easily applied for atom transfer between bands of 
even or odd parities, e.g. s ^ d or p ^ f and so on, 
and is therefore not suitable for preparing p-band atoms 
starting from the ground state. Surprisingly, the prop¬ 
erties of d-band bosons remain unexplored. Considering 
the recent experiment |28j . it is therefore of importance 
to work out the relevant physics of such systems. 

Here we analyse the zero temperature phase diagram of 
d-band bosons in an isotropic 2D square lattice. On the d 
bands we have three orbitals which are generally denoted 
da; 2 , dy 2 , and dxy However, due to the anharmonicity of 
the single site potential well, degeneracy is broken and 
the energy of the d^y orbital is somewhat higher than the 
other two orbitals which stay degenerate. The Mott in¬ 
sulator boundaries are determined within the Gutzwiller 
method. At low atomic densities, the d 2 ,y-orbitals are only 
weakly populated and we can therefor reduce the number 
of orbitals to two, i.e. we can describe the physics with 
an effective two-orbital BH model. The properties of the 
superfluid phase are analysed analytically by minimising 
the onsite energies. Like for the p-band superfluid phase, 
the order parameter is complex but due to the particu¬ 
lar shape of the d orbitals a new type of vortex lattice is 
formed with two vortex/anti-vortex pairs fixed to every 
site. At higher densities, the general structure of the vor¬ 
tex lattice persists but now also the da,y-orbitals becomes 
populated. 

One difference with respect to p-band BH models 
is the lack of a characteristic Z 2 -parity symmetry for 
bosonic atoms in the d band. The lack of this sym¬ 
metry stems from density-assisted orbital-changing col¬ 
lision terms. Nevertheless, there is another Z 2 -symmetry 
that is spontaneously broken in the superfluid phase. In 
the Mott insulating phase, we employ perturbation the¬ 
ory to derive an effective spin model, where it is seen 
that the new interaction terms appear in the form of an 
external field. In particular, the first insulating phase 
(the Mott with unit filling) can be described by an XYZ 
Heisenberg model in an external field. The model is, 
of course, still supporting the Z 2 -symmetry and we ex¬ 
pect a rich phase diagram within the insulating phase 
with the possibility of ferromagnetic phases with broken 
symmetry. It should be noted that the phase diagram 
of the 2D XYZ model in an external field is not fully 
known. We also present a brief discussion on how the so 


called Dzyaloshinsky-Moriya (DM) interactions (or anti¬ 
symmetric spin exchange) may appear when the lattice 
is no longer isotropic but the onsite degeneracy is kept. 
While not studied in detail, we point out that the insu¬ 
lating phases with higher fillings are also of great inter¬ 
est. Such phases will represent exotic magnetic states of 
higher spins or quasi-spins of SU{N) models. 


II. d-BAND BOSONS 

Throughout we will assume that all atoms reside on the 
d bands, i.e. the loading is perfect and that any time- 
scales are short compared to the decay times, driven by 
intra-atomic collisions, of these meta-stable atoms. The 
second condition implies that relaxation occurs on a time 
shorter than a few ps which is the experimentally mea¬ 
sured life-time of d-band atoms [55] . Relaxation typically 
happens during few tunneling times |24j , which allows co¬ 
herence on these excited bands to form before substantial 
dissipation/decoherence sets in. This has in particular 
been experimentally demonstrated for bosons loaded in 
the p and the / |57j bands. 

Before discussing the physics on the d bands, let us first 
recapitulate the main results for the system of bosons in 
the p bands of 2D square lattices. On the isotropic lat¬ 
tice, p-band bosons support two degenerate orbitals and 
a very interesting phase diagram. As for the regular s- 
band BH model Pi PH], the interplay between repulsive 
interaction and kinetic hopping leads to Mott insulating 
phases with integer fillings no = 0, 1, 2,... and condensed 
superfluid phases IHIIE]. The orbital degree of freedom 
implies, however, that both the insulating and the su¬ 
perfluid phases are non-trivial. In the condensate, the 
superfluid order parameter il){x,y) is a complex function 
representing a condensate with a vortex pinned to every 
site. The properties of the tunneling of p-band atoms 
in 2D further implies that vortices in neighboring sites 
have angular momentum quantised in opposite directions 
(i.e. winding numbers are ±1) |14H16j . Thus, there is 
a ‘checkerboard’ lattice of vortices and the two possible 
configurations underline the discrete Z 2 -parity symme¬ 
try of the model, as well as the spontaneous breaking of 
time-reversal symmetry. The orbital degree of freedom 
also makes the insulating phases very rich. The no = 1 
Mott phase with one atom per site can be mapped onto 
an effective Heisenberg XYZ model [T5] , whose phase di¬ 
agram in ID is qualitatively known [30] but not known in 
2D. In the spin language of the Heisenberg XYZ model, 
the breaking of the Z 2 -parity symmetry represents for 
example an Ising-like transition. 

The goal of the present paper is to demonstrate that 
also the experimentally relevant d-band bosonic model 
hosts very interesting physics with some notable differ¬ 
ences from the physics in the p band. A first guess would 
be that on the d bands, one again encounters a checker¬ 
board of quantised vortices but with winding numbers ±2 
instead. This would be the direct generalisation of the re- 
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suits appearing on the p bands. To construct such states 
all three orbitals need to be populated, but as we will see 
in the following only two of them are actually populated 
due to broken degeneracy. As a result we find a differ¬ 
ent type of vortex lattice. However, before discussing the 
possible phases on the d band, we give the corresponding 
many-body Hamiltonian in the next Section, and outline 
how we will analyse it within the Gutzwiller mean-field 
approach. 

A. Single particle properties 



Figure 1: The three d-bands; E^ 2 {qx,qy), Ey 2 {qx,qy), and 
Exy{qx,qy) for the potential amplitude V — 20Er. The band¬ 
width gives an estimate of the tunneling strength and the 
curvature the sign of the tunneling coefficient 

The single particle Hamiltonian for an atom with mass 
m in an isotropic 2D square optical lattice is ^ 

+ dl)+y [sin^ikx) + sm^{ky)] . (1) 

For convenience we introduce dimensionless parameters 
by letting the recoil energy = h?k^/2m set the en¬ 
ergy scale and the inverse wave number k~^ the charac¬ 
teristic length scale. Thus, in dimensionless variables, 

Hsp = ^ = -dl - dy +V [sin^(a;) -l-sin^(y)], with 

V = V/E-c- The spectrum E^{qx, Qy), consisting of bands 
separated by forbidden gaps, is characterised by three 
quantum numbers, a discrete band index v and two con¬ 
tinuous quasi momenta q^. and qy. In scaled dimensions, 
the first Brillouin zone is defined by qx, qy G [—1,-|-1). 
Increasing the potential amplitude V implies that the 
band gap increases while the widths of the energy bands 
shrink. The ‘flatness’ of the bands determines the mobil¬ 
ity of the atom. Consequently, as V increases the mobil¬ 
ity is reduced and the bands become flatter. In the limit 
of very deep lattice each site mimics a two-dimensional 
harmonic oscillator and hence the isotropic 2D case is 
characterised by a single s band, two degenerate p bands, 
three degenerate d bands and so on for the higher excited 
bands. The similarity with an isotropic 2D harmonic os¬ 
cillator can also be seen from the quantum numbers of 
the onsite orbitals; if \nx, ny) represents the eigenstate of 


the oscillator with quantum numbers Ux and Uy then the 
ground-state is the ‘vacuum’ |0, 0) (s-orbital state), the 
first excited states are |0,1) and |1,0) (p-orbital states), 
and the second excited states are |0, 2), |2,0), and |1,1) 
(d-orbital states). For finite, but large V, the separa¬ 
tions between the different set of bands (s, p, d ,...) may 
remain much larger than the band widths (at least for 
the lower bands). However, due to the anharmonicity of 
the potential wells the bands are no longer equidistant 
from each other, and even the three d bands may split. 
As an example, in Fig. we display the three d bands 
{Ex 2 ,Ey 2 , and Exy) for a potential of amplitude V = 20. 
As we can see, the dxy band has a higher energy than the 
other two. 



Figure 2: The three d-band Wannier orbitals of an isotropic 
2D square lattice. Here, blue/red colour represents nega¬ 
tive/positive values. Thus, the dx 2 and dy 2 orbitals have each 
two nodes in either the x- or the y-direction respectively, while 
the dxy-orhit has a single node in both directions. 

The eigenstates of the single particle Hamiltonian Hsp 
are the Bloch states with u = (rix,ny) the band 

index, q = {qx,<ly) the quasimomentum and x = {x,y). 
Taking the modified (i.e. restrict the integration to the 
first Brillouin zone) Fourier transform of these, one ob¬ 
tains the Wannier functions w,yi{x.) where i = {ix,iy) 
with ix, iy G JV labelling the site. Contrary to the de¬ 
localised Bloch functions, the Wannier functions are lo¬ 
calised around the site i and can therefore be normalised 
in the usual way. In addition, although they are not the 
eigenfunctions of the single-particle problem, Wannier 
functions are orthogonal and therefore provide an alter¬ 
native basis for describing such systems. In the infinitely 
deep lattice limit, they become simply the harmonic oscil¬ 
lator {x.\nx,ny) eigenfunctions localised at each site. For 
large but finite V, the harmonic oscillator eigenfunctions 
are not exact, but they still preserve the general struc¬ 
ture of That is, for the s band, where v = (0,0), 

Wsi(x) is approximately Caussian. For the p bands, of 
V = (1,0) and (0,1), the Wannier functions Wp^i{x) and 
Wpyi{'x.) are approximately Caussian in one direction and 
have a node in the transverse direction (the subscripts x 
and y tell the direction of the node). In the d bands, 
v = (2,0), (1,1), and (0,2), the three Wannier func¬ 
tions are shown in Fig. These are the three atomic 
d-orbital states, and the notation used is borrowed from 
atomic/chemical physics: The d 2 , 2 -orbital state has two 
nodes in the x direction and no node in the y direction. 
The opposite is true for the dy 2 orbital, while the dxy or- 
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bital has a single node in both directions. Furthermore, 
since the square lattice is separable, Wannier functions 
in 2D can be constructed as a product of ID Wannier 
functions. Explicit expressions for the d-band Wannier 
functions are: 


Wd^ 2 i{^) = Wsi{y)wdi{x), 


Wd^2ii^) = Wdi{y)wsi{x), (2) 


Wd^yii^) = Wpi{x)wpi{y), 

where we have introduced the subscripts s, p, d to label 
the Wannier functions of the different bands in ID. 

Naturally, the shape of the Wannier functions play a 
crucial role in the dynamics of the system: On a single¬ 
particle level, since the da orbital {a = y^) is much 

wider in the direction of the two nodes, tunneling pro¬ 
cesses parallel to this direction happen with considerably 
larger amplitude than in the perpendicular one. Thus, 
a (i 3 , 2 -orbital atom is more mobile in the x than in the 
y direction. The d^y orbital, however, is equally mobile 
in both directions. This direction-dependent propaga¬ 
tion results in an anisotropy in the problem, that is also 
present in the p-band system [I4HI6] (although absent on 
the s-band). This is also reflected as an anisotropy in the 
single-particle spectrum. As illustrated in Fig.j^ if 3,2 has 
a much larger curvature (inverse effective mass) in the x 
direction than in the y direction. Moreover, the sign of 
the curvature is also of importance since it determines 
whether the effective mass is positive (particle-like) or 
negative (hole-like), and as we will see below, it also sets 
the sign of the tunneling coefficient in the extended BH 
model. Furthermore, on a many-body level the ampli¬ 
tude of scattering processes is determined by overlaps of 
products of Wannier functions, and therefore their shape 
strongly influences the possible interaction processes as 
well as their strengths. 


B. Many-body Hamiltonian 


The Hamiltonian of the full many-body problem de¬ 
scribing contact interactions between the neutral atoms 
is given by 


^BH= / dx 


^^(x)ffsp^(x) -h^^^(x)^1'(x)^(x)^(x) 


(3) 

where the coupling constant Uq = Airti^a/m, with m the 
atomic mass and a the s-wave scattering length. The field 
operators ^ (x) and (x) annihilate and create a parti¬ 
cle at the position x and obey the bosonic commutation 
relations [4'(x), Tf (x )] = (5(x — x ). 


In order to constrain the atoms to the d bands (i.e., 
impose the single-band approximation)^ we expand the 


field operators in the corresponding Wannier functions. 


^(x) = ^ d2,2-,Wd^2-M) + dy2iWd^2i{^) + dxyiWd^yi{yd) 

. (4) 

with the sum running over all the lattice sites, 
and the bosonic creation/annihilation operators obey 

d/ 3 i, djjj = dij, where j3 = {x^, y^,xy} and all remaining 
commutators vanish. Following the usual prescription, 
^(x) and its Hermitian conjugate 'ii'f(x) are inserted 
into Eq. ([^. We then impose the tight-binding approx¬ 
imation such that the tunneling is restricted to nearest- 
neighbour sites and interactions to only dominant on¬ 
site terms. This yields the final form of the many-body 
Hamiltonian describing bosons in the d band. For later 
convenience, we separate it into two different parts, one 
containing all the processes that involve the dxy orbital, 
and one part that only include contributions of the d^^ 
and dy 2 orbitals: 


ffsH = idp + idd. (5) 

Furthermore, the two separate parts can be split up into 
dip — Hqp -f Hpkin T .f^pden T Hpc “f llpo^ 

( 6 ) 

Hd ~ Hod T .^dkin “f .^dden T Hdc T Udo- 

The free part of the Hamiltonian accounts for the onsite 
energies Ep of different bands together with the chemical 
potential p and is given by 

Hod = 'y ', {dd/3 tdjflfSi, (7) 

/3 


with the particle number operators n^i = d^-^dpi. The 
kinetic energy contributions are 


i7, 


pkin 


= -fP 'S^ 3} d ■ 
(ij> 


77dkin = — 

(ij).. 


( 8 ) 


with y2i{ij) the sum over nearest neighbours in the di¬ 
rection tr = {x,y}, y2{ij) the sum over nearest neigh¬ 
bours in all directions and a = {x^, j/^} the label of the 
orbital state. The tunneling coefficients satisfy 

\q\ > |tp| > \ti\, (9) 

where t^ and are the tunneling coefficients in the di¬ 
rections parallel and perpendicular to the two nodes of 
the orbital state. Furthermore tjj*, > 0 and fP < 0. 
Explicit expression of the overlap integrals are given in 
Appendix 
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The interacting part of the Hamiltonian Q is char¬ 
acterised by various processes. It contains the density- 
density interactions 


U, 


2 ^xyii'^xyi ~^‘^Upx (flx^i'^xyi~^^y^i^xyi^ 


I^pden — ^ ^ 
i 

.ffdden = ^ ^ [f^x^ii'^x'^i 1) T f^y^iif^y^i 1)] 


-{-‘ZUxy'^x^i^y^i f : 


the orbital changing interactions 


( 10 ) 


Uy.. 




ft dt 


-?t dt 


+ ( + <^y 2 ;dy 2 idxyida;yi j , 


+ (^xyi^xy-Ax^i^V^i + d)j2idL;da:yida:i/i 


Hdc — ^ + d]i2-,d]^2idx^idx'^i 

. . 

and finally, the density assisted orbital changing collisions 


Hpo — ‘2Upxy ^ jtla:yidy2j -\- d^2[^xyi^x^i^ (I^) 


the Hamiltonian is invariant under the transformation 

d’x'^i t dy2 j/ , 

dy2i —>• dx^i', (15) 


dxyi t d: 


xyi' ■ 


This swaps both the dx^ and dy 2 orbitals as well as the 
lattice indices ix and iy. Note that the interchange of 
indices is necessary due to the tunneling anisotropy of 
the dx 2 - and (iy 2 -orbital atoms ^ t^^). In the next 
section we will see that this symmetry can indeed be 
spontaneously broken to give interesting phases. We 
should also mention that in addition there are two more 
Z 2 -symmetries. Since the Hamiltonian only consists of 
quadratic or quartic terms, it is invariant under a sign- 
change of all operators. This parity symmetry can be 
further decided into two discrete symmetries, one which 
change the signs of the dx 2 and dy 2 orbitals, leaving the 
dxy orbital unaltered, and one with the reverse operation. 
These last two symmetries are trivial in the case when 


the dxy orbital is unpopulated as in Subsecs. HI B 1 and 

Imcl 


C. Mean-field approaches 

In the simplest mean-field analysis, the ground-state 
Id^o) of the many-body system can be expressed as the 
direct product |33] 


|4')mF = \^d2hipdyi,i’d^yilh 


(16) 


and 


dido — Uxxy ^ ) j^^a;2j {f^x^i T 'dy'^i) dy2\ (I^) 


T dy2\ {j^x‘^\ T d^y^i) dx2i 


(14) 


Other interaction terms vanish due to the symmetries 
of the orbitals. Like for the tunnelling coefficients, the 
various coupling constants can again be found in the Ap¬ 
pendix 

As pointed out above, the system of bosons in the 
p band supports a Z 2 -symmetry which can be sponta¬ 
neously broken to give a vortex-checkerboard condensed 
phase or a spin-flop phase in the superfluid or insulating 
phase respectively. This symmetry arises in that case be¬ 
cause p-orbital bosons in different orbital states scatter in 
pairs and therefore the number of particles in each orbital 
state is preserved modulo 2. On the d band, this con¬ 
served quantity is broken by the density-assisted inter¬ 
actions (12). There exists, however, another non-trivial 
discrete Z 2 -symmetry in the many-body Hamiltonian ([^. 
With the notation i = (ixfly), we let i' = (iyflx), and 


of coherent states; ^<i„ 2 i: V'dxyi)i = 

' 0 da,i|V'd,,. 2 ij and l/’da:yi)i “ 

'tpd^y-Md^ 2 h l/’<i^ 2 i> V'<ixyi)i- The classical energy functional 

is given by Acl 'l/'d, 2 i. V'd^ 2 ij'0dxyi = MF(^'|.ffBH|4')MF 
and the mean-field solution |4 ')mf is obtained from 
self-consistently minimising E^i V'd.yi 

with respect to the (complex) amplitudes 
and ipd^ i and the particle conservation constrain 

EiEol^doiP + Eil'^dxyip = iVtot, where Ntot is the 
total number of particles. Note, in particular, that 
since the Hamiltonian Hbh is normally ordered [34j . 
the energy expectation values are simply obtained by 
replacing operators dai, dxyi, dj^i, and dl^- by the cor¬ 
responding coherent state amplitudes V'dai; V'diyij V’d i) 
and respectively. This crude mean-field method is 
capable of giving a qualitative correct picture deep in 
the superfluid phase with a large or moderate number of 
particles per site, as will be utilised in Subsection |HIB| 
However, once the interaction is increased relative to the 
tunneling, the system becomes insulating and the above 
approach fails to predict such a phase: the method calls 
for improvement. 
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As a next step to advance on this coherent-state ansatz, 


(161 we consider intra-site particle fluctuations, but still 
without intersite correlations 


|'I')Gutz = (g)|<^i>i, (17) 


where every single site state is 

~ ^ ^ ^nl.2n^2n^y\'^x^ jn-xy)! (18) 

71^2 I'^xy 


with \nx'^,ny 2 ,nxy)i a Fock state with da; 2 -orbital 
atoms in site i and so on. Note that normalisation im¬ 


plies X) 


Tlx ,?1 x 


^Tlx 'hlyTLxy 


= 1 for every site i. The 


corresponding semi-classical energy functional is similar 


to the above, Ecutz 




= Gutz('l'|FTBH|^')Gutz, 


and the amplitudes c„^ 2 n,, 2 nxH obtained again via 
self-consistently minimisation of the energy functional. 
This procedure, called the Gutzwiller mean-field method, 
is particular for its ability of identifying the existence of 
Mott insulating phases HISS]. However, although the 
factorisation 0 becomes a more accurate approxima¬ 
tion in higher dimensions, it already provides a good es¬ 
timate of the insulating-superfluid phase transition in 


2D I3S]. 


III. PHASE DIAGRAM 


In this section we characterise different phases that can 
be expected for bosonic atoms in the d bands. To quan¬ 
titatively determine the phase boundaries of the insulat¬ 
ing phases we use the Gutzwiller method discussed in the 
previous section. For the p-band system El, the prop¬ 
erties of the superfluid phase can be understood from 
analysing the single site problem on a mean-field level 
and then considering how the presence non-zero tunnel¬ 
ing terms establish global long range order. Applying the 
same recipe for the d-band system fails for low atomic 
densities; it suggests that all three orbitals are populated 
while from the Gutzwiller analysis it follows that almost 
all population is distributed in the ^ 3,2 and dyZ orbitals. 
For this reason, when identifying the superfluid state in 
cases where the number of atoms per site ng ^ 5 with 
the coherent state mean-field approach, it is reasonable 
to constrain the analysis to account for only two of the 
orbitals. At higher densities however, the analysis must 
include all three orbitals. 

Since mean-field methods are less reliable for probing 
the physics in Mott insulator phases, this regime will be 
studied here in terms of an effective spin model derived to 
describe the limit where the magnitude of any tunneling 
coefficient |t| <C \U\, where U is the typical interaction 
strength. 


A. Superfluid-insulator boundaries 


The boundaries separating the Mott insulators from 
the condensed superfluid phase are determined, as men¬ 
tioned above, by the Gutzwiller ansatz wave function 
( |17[ ). The condensate order parameters for the three or¬ 
bitals are given by ifdf, = Gutz(^'Md ,3 |^')Gutz- Within the 
Gutzwiller method, the insulating phase is described by 
vanishing order parameters, while it is non-zero in the 
condensed superfluid phase. 


Before discussing the insulating-superfluid transition 
we make a brief detour on other possible phases absent 
in our case due to the imposed approximations. When 
the lattice amplitude is not too large, the larger width 
of the higher band Wannier functions may lead to non- 
negligible contributions from interaction between neigh¬ 
bouring sites. The effective model then includes nearest 
neighbour interaction terms, which if large enough (com¬ 
pared to the onsite interaction terms), may lead to su¬ 
persolid phases m- The existence of such a phase have 
been demonstrated for bosonic atoms in the p band |31] 
and occurs when the strength of interactions on neigh¬ 
bouring sites is as large as one fourth of the onsite inter¬ 
action strength. For 30 < H < 70 as considered in this 
work, however, the strength of interactions on neighbour¬ 
ing sites is always less than 1 % and such terms can be 
safely neglected. Nevertheless, it could be interesting to 
analyse the possibility of using the d rather than the p 
bands for realising such novel phases. As one would be 
required to relax the tight-binding, and most likely also 
the single-band approximation [35], this study goes be¬ 
yond the scope of the present work. 


Let us now turn to the Gutzwiller predictions for 
the transitions. The self-consistent minimisation of the 


^rton.orix 


is based 


Gutzwiller energy functional F^Gutz 

on the use of the Nelder-Mead simplex algorithm [39j . 
Since we are interested in the phases for the system with 
different densities, the chemical potential y (independent 
of the orbital flavour) is added to the Hamiltonian, see 
Eq. Q. For the numerics we truncate the number of 
Fock states in the sum (18) to four for every orbital, i.e. 
up to three atoms for each orbital state and nine atoms in 
total per site. This gives an onsite Hilbert space dimen¬ 
sion of 64, and we can accurately capture the first four 
Mott lobes. Going to higher insulating states would re¬ 
quire Fock states of larger particle number, and thereby a 
rapid slow down of the numerics. Furthermore, we scale 
Eq. ([^ with the interaction strength U. 


The results of the numerical Gutzwiller calculation are 
presented in Fig. showing the superfluid order param¬ 
eter in the yt-plane. In this range of the system’s pa¬ 
rameters, the occupation of the d^y orbital (fixy) ~ 0 
everywhere. Throughout, proper Wannier functions are 
used in the computation of the overlap integrals ( |Al[ ) pre¬ 
sented in the Appendix]^ and we consider the potential 
amplitudes 30 < E < 70. These large potential ampli¬ 
tudes (up to 70 recoil energies) are required to ensure 
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Figure 3: The order parameter ip = |(dj. 2 i)| = |(dj, 2 i)| 
(The corresponding order parameter for the djjy-orbital is 
approximately zero except for t/U < 10“^ where it is, 
however, at least one order of magnitude smaller than 
As is seen, the chemical potential is varied such that the 
first four Mott insulators are illustrated (dark blue regions 
representing a vanishing superfluid order parameter). For 
this plot, the relative strengths between the interaction 
terms has been taken as {Uxy, Up, Upx, Upxy, Uxxy}/U = 
{0.17, 0.9, 0.3, 0.04, —0.03} corresponding to a lattice with 
amplitude V = 40. With this lattice depth, 75^.2 = Eyi = 0 
and Exy = 1.6[/, and the relative tunnelings , |t^|}/t“ = 
{0.002, 0.007}. 


the validity of the single-band and tight-binding approx¬ 
imations in the d band, that accommodates very broad 
Wannier functions. The order parameters for the 
and dy 2 orbital are identical and for simplicity we call it 
Ip = |((ia;2i)| = \ {dy2p)\. Note that the effective Gutzwiller 
tunneling t used in the figure is twice the sum of the s- 
and d-band tunnelings. More precisely, out of the four 
nearest-neighbours in the 2D system, tunneling processes 
in the d band occur with (ID) d-band tunneling strength 
to two of the neighbouring sites, and with (ID) s-band 
tunneling strength on the remaining two neighbouring 
sites. The general structure of the Mott lobes are similar 
to those found for the regular one orbital BH model us¬ 
ing mean-field methods. That is, the extent of the higher 
lobes fall off as roughly Hq^ for the filling factor no [5S] . 

For higher densities no, the typical interaction energy 
becomes considerably larger than the energy gap separat¬ 
ing the Exy from the and Ey 2 bands, and therefore 
the occupation of the d^y orbital is expected to increase. 
This regime can be explored numerically but requires 
a higher truncation in the number of particles, which 
makes computations much more costly. To circumvent 
this issue, we studied the system at fixed chemical po¬ 
tential and for various values of t/U and indeed, 
the results show that with the the same parameters as in 
Fig.[^ but with for example {^, t}/U = {5, 0.5} the pop¬ 
ulation of the three orbitals are (( 77 . 3 , 2 ), (hj, 2 ), (hjjy)) « 
(3.3, 3.3, 1.7). 


B. Onsite superfluid states 

The Gutzwiller mean-field method of the previous Sub¬ 
section provided us with an estimate of the Mott lobes, 
but the orbital dependence within the superfluid phase 
was left implicit. We will approach this problem next 
using of the coherent state mean-field ansatz for clarity. 


1. Moderate atomic densities 


Generalising the results from the p bands [mile] we 
would expect a doubly quantised vortex on every site for 
the d-band system. Such a state takes the form 


^'(x) = v^/2 


(x) - (x) ± V2iwd^y (x) 


<t4> 


^ = y/N,l2 


1 

-1 

±-\/ 2 * 


(19) 

where Ng is the number of atoms at the site and is fixed 
a priori. However, according to the previous Subsection, 
the occurrence of a state as that of Eq. (191 is precluded 


for lower atomic densities by the negligible occupation 
of the dxy orbital. Moreover, from only Wd 2 (x) and 
Wd 2 (x) it is not possible to construct a doubly quan¬ 
tised vortex state, which leads to the conclusion that the 
superfluid order parameter is a state of a different type. 

Following the results of the previous Subsection and 
focusing first on the low-density case, we approximate 
ipdxy = Oj rewrite the single-site order parameter as 


XiJ = 


ipd 


= y/K 


cos 9 
sin 9 


( 20 ) 


where 0 < 9, < tt and 0 < (p < 27r. Even though this is 
valid for systems with onsite particle number limited to 
10 , we believe this gives a good picture of the superfluid 
state. Indeed, as we show in the next Subsection, a sim¬ 
ilar state is found when computation includes the third 


orbital state. With the parametrisation of Eq. (20), the 
energy functional becomes 


^ = ^cos4d+^sin^20Q+cosV) 


( 21 ) 


+Uxxy sin 29 cos p, 

which after minimisation yields the fixed point {9o,po) 
9o = 7 r/ 4 , (po = arccos {-Uxxy/Uxy), (22) 

where dgEi,i[9o, Po] = Sipp'd [do, Po] = 0. This solution. 




(x) = U^°Wd^A:^) + Wd^^{yi) , (23) 

















is depicted in Fig. (a) and (b), for the parameters cor¬ 
responding to a lattice with V = 40. The first plot (a) 
shows the atomic onsite density |'Pvor(x)P and the second 
(b) the phase (p of the order parameter. As is evident, 
the onsite condensate hosts two vortex/anti-vortex pairs. 
Each vortex has a phase winding of ±27r around each sin¬ 
gular point. As argued above, while the naive harmonic 
approximation on the d bands suggests a single multiply 
quantised vortex in each site, they do not occur here due 
to broken degeneracy of the orbitals. It has also been 
known that doubly quantised vortices should be energet¬ 
ically unstable in a harmonically trapped Bose-Einstein 
condensate m- 

Throughout, proper Wannier functio ns a re used in the 
computation of the overlap integrals presented in 

the Appendix and we consider the potential ampli¬ 
tudes 30 < E < 70. These large potential amplitudes (up 
to 70 recoil energies) are required to ensure the validity of 
the single-band and tight-binding approximations in the 
d band, that accommodates very broad Wannier func¬ 
tions. Typical ratios between the interaction strengths 
for E = 40 are presented in the caption of Fig. 

The condensates at neighbouring sites have the same 
ordering of the vortex pairs since t± > 0, which sup¬ 
ports a “0” or “27r” phase locking between consecutive 
sites. Note that this is always the solution in case of 
occupation of only the d^^ and dy 2 orbitals, and that it 
does not depends on U. In addition, the configurations 
defined by the vortices orientations spontaneously 
break the Z 2 -symmetry defined in Eq. (151 as well 


as time-reversal symmetry. For the solution given in 
Fig. (a) and (b), and starting from the vortex in the 
upper left corner, the phase winds clockwise/counter- 
clockwise/clockwise/counter-clockwise. The alter¬ 
native configuration with broken Z 2 -symmetry 
has instead counter-clockwise/clockwise/counter- 
clockwise/clockwise winding. 

Now the natural question is what are the properties 
of the phase with restored Z 2 -symmetry? To answer 
this, we first notice that according to (22), the shape of 
the two-vortex pairs depend on the ratio R = Uxxy/Uxy 
When i? —>■ 1 the vortex/anti-vortex approach each other 
and annihilates at i? = 1. For R > 1, the order parame¬ 
ter reads instead 


^'sol(x) = Y ^ 


Wd 2 (x) ± Wd (x) 


(24) 


The atomic density and the order parameter phase are 
shown in Fig. (c) and (d). The density vanishes iden¬ 
tically around a circle, where also the phase makes a tt- 
jump from (po = 0 to ipo = tt. Such a behaviour describes 
a dark soliton, immobile and confined to the lattice site. 

If the ratio R could be externally controlled, then the 
system could be driven through a phase transition from 
a ‘soliton’-superfluid into a symmetry broken ‘vortex’- 
superfluid. In addition, since the derivative of the classi¬ 
cal energy is discontinuous at the critical point Rc = 1, 
it suggests that the transition is of the first order. This 



Figure 4: The onsite density |T(x)|^/As (a) and (c) and the 
corresponding phases Arg[4'(x)] (b) and (d). Here we as¬ 
sumed that the atomic density is rather low such that we can 
neglect any population of the d^y orbital. In the upper two 
plots we consider a lattice with an amplitude E = 40. It 
is seen in (b) that the phase of the order parameter winds 
47r at the four points where the density vanish. This reflects 
the presence of four vortices - two vortex/anti-vortex pairs. 
In the lower two plots we use the same parameters but put 
Uxxy = —2Uxy in order to reach the regime where the state 
qualitatively changes. Here a dark soliton is separating the 
central peak from the surrounding circle. 


is not so surprising, since the states separated by the 
critical point are topologically different. How to control 
R experimentally may be non-trivial, but nevertheless, 
the present analysis sheds some light on the underlying 
physics of d-band bosons. 


We end this Subsection by noticing that on the p band 
in 3D the orbitals are also three-fold degenerate. In this 
case an artificial instability appears whenever the har¬ 
monic approximation is considered, i.e. by replacing the 
Wannier functions by the corresponding harmonic oscil¬ 
lator eigenfunctions, which implies that it can become 
favourable to populate the orbitals unequally m El- 
On the d band, the depopulation of the dxy orbital does 
not derive from such a spurious effect, but instead, it 
appears together with the self-consistent minimisation of 
the full (lattice) energy functional, that includes both 
the interacting and tunneling processes simultaneously. 
This means, furthermore, that the combined (indepen¬ 
dent) knowledge of onsite properties with additional con¬ 
siderations accounting for effects of the tunneling - as one 
usually proceeds for studying the system in the p band, 
is not enough to determine the physics in the d band. 
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Figure 5; The same as Fig. |4]but deep in the superfluid phase 
(no ^ 1) where all three orbitals are considerably populated. 
The difference between the upper and lower plots is the sign 
of the order parameter (the two states have equal en¬ 

ergy), which reflects the states in two neighbouring sites in 
the lattice. Despite the fact that the d^y orbital is largely 
populated in this case, the general structure of the superfluid 
state shown in Fig.|^(a) and (b) survives, i.e. the onsite order 
parameter hosts two vortex/anti-vortex pairs. However, the 
state gets distorted with two of the vortices at the very edge 
of the distribution. As for Fig. the potential amplitude 
V = 40. 


2. High atomic densities 


In this Subsection we consider a larger particle number 
uq ^ 1 such that all three orbitals are populated. Re¬ 
laxing the assumption ipd^ = 0 adds a new variable to 
be accounted for in the minimisation of the energy func¬ 
tional. By using the parametrisation of the three-orbital 
order parameter in spherical coordinates. 


= -s/^ COS 9 cos 

'ipd ^2 = cos Osin (25) 

V'd.j, = v^sind. 


the normalization constraint is automatically taken care 
of and by further choosing the overall phase to be zero we 
are left with four angle parameters; Ec\[9, (j), ipj-d]. Since 
in the general case we do not find analytical solutions 
for the minimisation (fix point) the system is studied 
numerically (again utilising the Nelder-Mead simplex al¬ 
gorithm) . 

When occupation of the d^jy-orbital increases, it is in 
principle possible to form doubly excited vortices in each 
of the lattice sites. However, as argued above, it is 
questionable whether such state could be energetically 
favourable. Another aspect arising with all three or¬ 
bitals populated is that the other two Z 2 -symmetries 
discussed after Eq. (15) become relevant as will be 


demonstrated. These two symmetries are reflected in 
the fact that the energy Ed is invariant under ei¬ 
ther {itd^2,i’d^2,'4’d^y) ^ (.-^ 2 ,-^d^ 2 ,^d^y) and/or 
ii’d^ 2 ,^d^ 2 ,^d^y) ^ (.i’d^ 2 ,'^d^ 2 ,-'^dj) (we note that 
these two are equivalent via a gauge transformation of 
the overall phase). 


The numerically obtained ground state is displayed in 
Fig. [5 (a) and (b). Here, the population of the three 
orbitals dx 2 , dy 2 , and dxy are 0.36, 0.33, and 0.31 respec¬ 
tively, and hence none of the orbital states dominates the 
others. The effect of the non-zero population of the d^y 
orbital is to squeeze and rotate (in this example counter¬ 
clockwise) the atomic distribution without destroying the 
two vortex/anti-vortex pairs; the larger population in the 
dxy orbital the more rotated and squeezed is the atomic 
distribution. Thus, as is clear from the deformations 
in Fig. vortices of (for example) positive winding are 
closer to the centre of the site where the density is higher. 
This suggests that as interactions become stronger and 
the different orbitals become effectively more degenerate, 
the squeezing is also increased and finally two of the vor¬ 
tices appear were the density is vanishingly small. In this 
case the onsite states start to approach doubly quantised 
vortex states. This is an interesting observation, since as 
we have pointed out earlier such a vortex is typically not 
stable. However, one should bare in mind that here we 
have two singly excited vortices coming infinitely close 
to each other. Furthermore, this continuous deformation 
does not imply a violation of conserved angular momen¬ 
tum; the state only appears as a doubly excited state, 
and in the full lattice model the doubly excited vortices 
would alternate between neighbouring sites such that the 
total angular momentum vanishes. 


To understand the full superfluid state that extends 
over the whole lattice we recall that the tunneling am¬ 
plitude H on the p band comes with a negative sign. 
This implies that ipd^^y changes sign between neighbour¬ 
ing sites while 'ipd ^2 4’dy2 do not. This is illustrated 
in Fig. (c) and (d). The reversed sign of il^d^y has 
the effeCT of rotating the distribution in the opposite 
direction (i.e. here clockwise), but the winding num¬ 
bers are unaffected. Now recall that the Z 2 -symmetry 
of Eq. (15) corresponds to reversal of the phase wind¬ 
ing. Taking all into account we have that in the symme¬ 
try broken phase the onsite rotation direction alternates 
between neighbouring sites (breaking of the parity sym¬ 
metry i^d^ 2 ,^d^ 2 ,'^d^y) ^ {^d^ 2 ^^d^ 2 ^-^d^y)) and fur¬ 
thermore the phase winding is determined from breaking 
of the symmetry (15). In addition, we have numerically 
checked that in the same way as in the previous subsec¬ 
tion, the state in the phase with restored symmetry is 
characterised by a dark soliton. 
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C. Insulating phases for filling no = 1 

Deep inside the Mott insulating phases the relevant 
physics of the many-body system can be studied in terms 
of an effective (pseudo) spin model derived from the per¬ 
turbative treatment of the tunneling processes relative 
to the interaction |42] . By mapping the problem onto a 
spin system, the effective dimension of the Hilbert space 
is greatly reduced and the relevant physics becomes more 
transparent and easily analysed. Moreover, even though 
the properties of many spin models are still unknown, 
known results of previous studies can be applied to the 
present case. 

Since the different Mott lobes are characterised by a 
definite integer filling no, the derivation of the effective 
spin model can be most naturally handled with the use 
of operators that project the eigenvalue problem onto or¬ 
thogonal particle number subspaces of the Hilbert space. 
As a consequence, the result of the perturbative treat¬ 
ment is not general to any insulating state, and thus we 
restrict the study here to the first Mott lobe (no = 1) 
which is most relevant for experimental studies. We pro¬ 
ceed by defining the P and Q operators, = P and 
= Q, that project, respectively, onto the space of 
states with unit filling, and the states with at least one 
doubly occupied site. For the lowest Mott insulator we 
may again impose the assumption that the cia;y-orbital 
state is negligibly populated, and that the only relevant 
degrees of freedom considered stem from the dj ;2 and dy 2 
orbitals. 

The eigenvalue problem can then be written as (for a 
more detailed derivation see the Appendix and Refs. [THl 


Hu = Hdden + Hdc + iFdo, cf. Eqs. ([^-([^. Due 
to the assumptions of the Mott phase, the resolvent 

1 / {qHvQ - e) can be expanded such that the effec¬ 


tive Hamiltonian contains contributions to second order 
in t/U. In addition, the tight-binding approximation al¬ 
lows proceeding with the analysis as a two-site problem. 
Therefore, we define the basis spanning the P subspace 
by 


Hv = { \x, x), \x, y), \y, x), \y, y)}, (27) 

where \a,l3) = corresponds to the state with a 

da-orbital atom in the site i and a d^-orbital atom in the 
neighbouring site j, a,/3 = {x, y}. In the same way, the 
relevant states in the basis of the subspace Q of doubly 
occupied sites is given by 


Hc = {|0,2x),|0,xy),|0,2y)}, 


(28) 


with |0,2a) = 2 and |0,a/3) = dijdjsjlO). 

Due to the possibility of transferring population be¬ 
tween the different orbital states via H^c and Tfdo, we 
notice, however, that the projection of the Hamiltonian 
onto the Hq subspace is non-diagonal in the basis of in¬ 
termediate states of the perturbation theory. A practical 
way of dealing with this situation is to notice that the 
contributions of the different processes can be obtained 
from computation of {Hq — E), where Hq = QHtjQ, 
with a subsequent inversion. In addition, since the en¬ 
ergy E ~ jU we have [Hq — E)~^ 


Hq ^. More 


explicitly, ordering the basis oi Hq according to (28), 


Hyiq — i — P HdkinQ 


1 


QHuQ - E 


QiJdkinP, (26) 


Hq = 


Uxx 

V2Uxy 

Uxy 

V2Uxy 

2Uxy 

V2Uxy 

Uxy 

V2Uxy 

Uyy 


(29) 


where Hdkin is given by Eq. 


and the interaction part 


and 


= A 


2UxyUyy U^y 

V2{U^y - U^yUyy) 
J2 


-2H 


xy 


y/2[Uly - U^yUyy) 
2U,,Uyy - 2Uly 
^{Uly - U^yU,,) 


-‘^U^y 

V^iUly - U,yU,,) 

‘2Ua:yU^x - Uly, 


(30) 


with 


A = {PU^.UyyU^y - 2U,,U^y - 2UyyU^yy 


(31) 


by computing the relevant matrix elements of Eq. (261. 
This yields (see Appendix for more details of the 
derivation) 


We determine the final form of the effective Hamiltonian 










11 


Hno = l — ■ 


- Ulp) n„in„j + 2|t“|2 {U^^Upp - U^p) n^,hp; 

a/3 (ij)^ 

-2t°t^Uapd^pidaidl^daj + 2t°t^ {UaaUpp - U^p) rfjjjdaidjjjd/3j 
+ (|i“P + Cta) {Ulp - UapUpp) flaidl-dcj + {U^p - UapU^a) 4i^^/3i^aj } 


(32) 


By further employing the Schwinger angular momentum 
representation [15] 


4+ = St + iSf = (33) 

5r = Sf - iSV = dl,,d,2, 


together with the constraint of unit filling ft; = h^ 2 i + 
hy 2 i = 1, Eq. (321 can be mapped into a spin-1/2 XYZ 
model with DMlnteractions |44l|45] and an external field 


Hxyz 


-E^ 

(ij> 


(l + 7)5r5/ + (l-7)^f^]' 


[AS?St + <5 {StS? + Sts?)] + TS?. 
(ij> i 

(34) 

Here J = with U"^ = 4:{UxxUyy — U^y) and 

A given by Eq. (31), the anisotropy parameter 7 = 
-4C/JJC/2, 


A = \tlW2UxyUyy-Uly) + \tl\^2UxyUxx-U, 
-{\tl\^ + \tl\^) {Uly+WxxUyy), 


^ 1 
xy) 


5 = {Uly-UxyUyy){\t%\^+t%ty) 

-{uly-UxyUxx){\tl? + t%ty), 

r = {uly-UxyUyy){\ti\‘^ + tity) 

+ {uly-UxyUxx){\ti? + t%ty). 


However, in the isotropic lattice as considered in this 
work, the DM interactions vanish due to the Z 2 - 
symmetry discussed in Eq. (15), which in the spin lan¬ 
guage read S? —>■ Sf, Sy —>■ —5/, and S? -P- —S?- This 
symmetry is broken in anisotropic lattices. In fact, there 
ty ^ ty and tt_ and in such case, since <5^0, the 

DM terms reappear in the effective spin model. 

Several interesting facts should be noted. First, that 
in the same way as for the system of bosons in the p 
band m, d-orbital bosons provide an alternative reali¬ 
sation of the XYZ Heisenberg model, albeit here with the 


presence of an external field even in the isotropic case. 
This is a consequence of the fact that the density assisted 


processes, Eq. (12), break the Z 2 -parity symmetry asso¬ 
ciated to orbital changing interactions of Eq. (11) which 
preserves the number of orbital atoms modulo two. This 
also explains the presence of the Sf term in the Hamil¬ 
tonian and it readily follows that in the absence of this 
field the parity symmetry is restored. Second, that in the 
same way as for the external field, the DM interaction re¬ 
sults from the density-assisted processes. Third, that in 
the limit of vanishing density-assisted interactions, the 
effective spin Hamiltonian becomes identical to the cor¬ 
responding one of the p-orbital system, of Ref. [15]. That 
this indeed should be the case can be understood from 
the fact that the anisotropy parameter 7 is a consequence 
of the orbital changing processes. However, a difference 
between the two models is that while the t^t^ < 0 for the 

p-orbital system, f^t^ > 0 in the d band, and therefore 
the effective spin model favours primarily ferromagnetic 
alignment at neighbouring sites in all the pseudo-spin 
components. 


To the best of our knowledge, the phase diagram of 
this model is not fully known in 2D even for the simplest 
case of the isotropic lattice, where <5 = 0. However, the 
underlying Z 2 -symmetry suggests the possibility of a rich 
phase diagram. We note that for T/J A/J (I 7 I < 1) 
the system is characterised by a highly magnetised state 
in the a:-spin component, while in the opposite limit it is 
ferromagnetic in the z-component of the spin (noting that 
A > 0). The later is a symmetry broken phase with mag¬ 
netisation in either positive or negative z-direction. For 
J ^ 7, A another possible symmetry broken ferromag¬ 
netic phase could appear with the spin in the a/y-plane or 
possibly a gapless floating phase which exists in the ID 
XYZ model. In ID the qualitative features of the XYZ 
phase diagram are known m- This has been studied 
in terms of p-band bosons [18] . but following the same 
procedure of that study to reduce the dimensionality to 
an effective ID model one would inevitably generate the 
DM interaction terms for the system in the d band. This, 
however, gives interesting possibilities since this model is 
known to host an interesting phase diagram with ferro¬ 
magnetic and Luttinger liquid phases [46] . 












12 


IV. CONCLUDING REMARKS 


Motivated by a recent experiment [55], we have stud¬ 
ied the zero temperature properties of bosonic atoms on 
the d bands of an isotropic square optical lattice. Due to 
the particular shapes of the onsite d orbitals, the phases 
characterising the ground-state of this model are differ¬ 
ent from those analysed in the past for systems on the 
p band MM- This is not surprising if we think of the 
orbitals in the harmonic approximation where the p or¬ 
bitals carry angular momentum components Z = ±1 and 
the d orbitals I = 0, ±2. The |Z| = 1 angular momentum 
of the p orbital is reflected in the singly excited vortex 
at each site in the superfluid phase. The higher angular 
momentum components of the d orbitals implies that the 
onsite state can take different forms than a state with sin¬ 
gle vortices. The direct generalisation of the superfluid 
phase on the d band would be a checkerboard lattice of 
a single |Z| = 2 vortex located at every site. However, 
we found instead that two vortex/anti-vortex pairs ap¬ 
pear on each site. The common wisdom is that multiple 
excited vortices are energetically unstable to form sev¬ 
eral singly excited vortices |40| , and one could argue that 
this explains why we should not And doubly excited vor¬ 
tices in this system. This simplified argument cannot be 
true here since the total angular momentum should be 
preserved, which is not the case for a state carrying two 
vortex/anti-vortex pairs where the total angular momen¬ 
tum vanishes. 


The inter-site phase locking of the different intra-site 
order parameters is determined by the signs of the dif¬ 
ferent tunneling coefficients. For the d band this implied 
that the corresponding single site vortex states are or¬ 
dered in the same way in all the sites. This is in con¬ 
trast with the checkerboard type (or anti-ferromagnetic) 
arrangement of vortices in the superfiuid phase of the 
system in the p band. This ferromagnetic phase is a 
symmetry broken phase with respect to the direction of 
phase windings of the vortices. The spontaneous break¬ 
ing of the symmetry (151 is also accompanied by breaking 
of time-reversal symmetry, which otherwise would imply 
the possibility of switching between the two degenerate 
vortex states of different orientation. In the symmetry 
preserved superfiuid state, a dark ring soliton is formed 
on every site. This phase may, however, not occur nat¬ 
urally in experiments since the strength of the density 
assisted interaction coefficient must become comparable 
to the coefficient of the orbital changing interaction. This 
said, it does not necessarily mean that this phase cannot 
be monitored but then by external driving like suggested 
in Ref. [T8l. 


Like for the superfluid phase, the orbital structure also 
makes the physics of the insulating phases very intrigu¬ 
ing. We focused on the lowest Ailing no = 1, where the 
Mott phase was shown to feature very interesting prop¬ 
erties. With a perturbative treatment, this system was 
mapped onto an XYZ Heisenberg spin model in an exter¬ 
nal held. The phase diagram of this model is not known 


in 2D, but by considering limiting cases we argued that 
the model should posses a very rich phase diagram. We 
further showed the appearance of DM interactions in the 
case where the lattice is no longer isotropic. In particular, 
by fine tuning the lattice is is possible to restore the de¬ 
generacy between the dx 2 and dy 2 orbitals and still break 
the lattice isotropy. In the harmonic approximation this 
accounts to fulfilling the condition Vxkl = Vy ki im ( in 
unsealed variables, where Vy and are the potential am¬ 
plitudes and wave numbers in the directions p = x, y). 
Thus, it is possible to study effects arising from the DM 
interactions, of which the most famous effect is perhaps 
that of spin canting where an anti-ferromagnetic state 
builds a finite magnetisation or where the magnetisation 
in a ferromagnetic state gets quenched [JS] . The presence 
of these terms imply the breaking of the Z 2 -symmetry 
(151 and as a result the transitions should turn into first 
order or of Berezinskii-Kosterlitz-Thouless type [l7] . 


In 3D, the no = I Mott insulator could be effectively 
described by an SU{3) pseudo-spin model, as was re¬ 
cently shown for the system in the p band [21|. T he den¬ 
sity assisted orbital chaining collision terms (12) present 
on the d band would, however, give some additional terms 
to the effective model in comparison to that emerging on 
the p band. Now, what could we expect from the physics 
on other insulating phases? On the ng = 2 insulating 
phase, the projected (two atoms/site) single site Hilbert 
space is spanned by the three states {|a;a;), \xy)^ \yy)}- 
In the Schwinger spin representation ( |33| ) with the con¬ 
straint hi = 2 one would obtain an effective spin-1 XYZ 
model. The integer spin implies the possibility of a gap¬ 
less Haldane phase to appear |48|. Going up into the 
Mott lobes with larger filling, the effective spin models 
would be characterised by higher spin, but as long as we 
remain in 2D, the (pseudo) spins in the two-orbital case 
would be generators of the SU{2) group. If occupation of 
the dxy orbital becomes non-negligible, then the pseudo¬ 
spins are the generators of the SU{3) group even in the 
2D lattice, with a corresponding spin model with DM 
interactions in all the components. 


Since this work presents the first theoretical study of 
the physics in the d band, we have limited the analysis 
to the isotropic 2D lattice. Naturally, other lattice ge¬ 
ometries and in other dimensions the physics is expected 
to change as well as if one considers fermionic atoms in¬ 
stead of bosonic ones. Another aspect left out here is the 
influence of a trapping potential. On the p band it was 
recently demonstrated that the presence of a harmonic 
trap is a most simple way to directly detect outcomes of 
the anisotropic tunneling on the excited bands [231 • Sim¬ 
ilar effects as those found for the confined p-band sys¬ 
tem would also appear on the d band, i.e. time-of-flight 
detection of the freely expanding atomic cloud or single¬ 
site addressing |50] would reveal information about the 
intrinsic anisotropic tunneling on the d bands. 
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Appendix A: Hamiltonian parameters 


Here we give the general expressions for the overlap 
integrals rendering the various parameters of the many- 
body Hamiltonian (§. Making use of the separability of 
the potential and the Wannier functions ([^, and after 
introducing the ID single particle Hamiltonian Hisp = 
—dl + V sin^(a;), the coefficients become 


t 


S _ 

± — 


dx Wsi+l {x)HispWsi (x), 


• To joi, ttj) via action of which contribute to 

the effective Hamiltonian with terms of the type 

-EE 2^ (Bl) 

(ij)^ “,/3 


• To via action of 

-EE 2^ (C^-C//3/3 - nai^jdcj. (B2) 

(ij><T a,/3 


tP = 


dx w*i+i {x)HlspWp^ (x ), 


To via action of dLdygj; 


ty = y dxWdi+l{x)HispWdi{x), 

^ ={^I d-xw%[x)^, 

Up = (^J dx\wpi[x)\'^ , 

Uxy= (^J dx w%{x)w'^,i{x)^ , 

Uxxy=(J dxw%{x)wsi{x^(J dxwdi{x)w%{x^, 
Upxy= (^j dx\Wp,{x)\'^Ws^{x)Wd^{x) 


Upx = (J dxwli{x)\wpi{x)\‘^(J' dx\wsi{x)f\wpi{x)\‘^. 

(Al) 

In the above expressions we have used the phase conven¬ 
tion that Wsi{x) and Wdi{x) are purely real while Wpi{x) 
is purely imaginary. 


Appendix B: Explicit computation of the effective 
spin model 


.tZtd 


- E E 2if - Ulp) d%d. 


(ij)<T “A 


ai^ay (B3) 


To |Ai/3j) via action of 

E E 2~^^a/3 d^idai^^jdaj. 


(ij)a a./3 


(B4) 


The states of the the type |ai,^j) are also connected via 
tunneling to the three intermediate states in the Hq sub¬ 
space; 


Kdl.daj\ai,l3j) = A:|0,aj^j) = 

(a“|| 0,2aj) + A“^|0,aj/3j) + 2/3j)) . 


(B5) 


He re, in addition to the conjugates of Eqs. (B2) 
and (B3), the other possible transitions are 


• To joi, /3j) via action of d^iday which results in the 
contribution 

“ E! 2 A {UaaUpp — U^p) flainpj. (B6) 


(ij>,T a./3 


To via action of d}p-Jipy, 


The final form of the effective Hamiltonian is obtained 


after computing the relevant matrix elements of Eq. (26). 


That is, we consider all the different transitions allowed 
for each state. 

The states of the type |Q;i,aj), i.e. the same type of 
orbital atom in the two neighbouring sites, are connected 
via tunneling to the three different intermediate states in 
the Hq subspace; 

k dj^jdailoi, a-^)y/2k\Q, 2aj) 


= V2 (if““|0, 2a;) + |0, aj/3j) + Add|o, 2/3j)), 

where we have introduced the shorthand notation K = 
Hq^. The possible transitions are 


- E E 2^ - Ulp) d^L-A.dp;. (B7) 

(ij><T «:/3 


Combining the matrix elements of Eqs. (B1)-(B7), we 
derive the effective Hamiltonian (321. 
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